%%%%%%%%%%%%%%%%%%%%%%%%%
% SIMULATE INDIVIDUALS %%
%%%%%%%%%%%%%%%%%%%%%%%%%
% note, when passing iX, forcing a given draw of states 
% for replication purposes

function simulation_output = solve_main(model_parameters,iX)

nInd = model_parameters.nInd;
CAPH = model_parameters.CAPH;
MU = model_parameters.MU;
rho = model_parameters.rho;
%r= model_parameters.r;
income_file = model_parameters.income_file;
gamma = model_parameters.gamma;
capHtmod = model_parameters.capHtmod;

T = size(MU,3);
Nstates = size(MU,1)-1;

% simulate histories
% start with the capHtmod distribution

if ~exist('iX','var')
    iX0 = randomDiscrete(capHtmod(ones(1,nInd),:)');
    iX = zeros(T,nInd);
    iX(1,:) = iX0;
        for t = 2:T
            capHt = CAPH(:,:,t-1);
            iX(t,:) = randomDiscrete(capHt(iX(t-1,:),:)');
            % iXt gives col of fm matrix
            % draw row from a 1-100 uniform
            % for each individual
        end
end

% table 2 (simulated)

table2 = [];
for k = 1:(T/5)
    l = 5*(k-1)+1;
    age = l + 24;
    if age==75
        h = T;
    else
    h = l+4;
    end
    delta = h-l+1;
    a = tabulate(reshape(iX(l:h,:),delta*nInd,1));
    a  = a(:,3);
    if length(a)<(Nstates+1)
        a = cat(1,a,0);
    end    
    table2 = cat(1,table2,a');
    if age==75
        break
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate german contracts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

maxreps = 100;
tolFP=0.005;
Pspot_N = PGR_prices_v2(MU,CAPH,rho,T,tolFP,maxreps);

%[Ppaid_D Ppaid_mean_D Switch_mean_D dead_mean_D Stay_D] = Paid_prices(Pspot_D,iX);

[Ppaid_N Ppaid_mean_N Switch_mean_N dead_mean_N Stay_N] = Paid_prices(Pspot_N,iX);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VI. prices paid in short term (full) contracts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = squeeze(MU);
t = repmat((1:T)',1,nInd);
Ppaid_spot = a(iX+(Nstates+1)*(t-1));
Ppaid_spot_mean = sum(Ppaid_spot,2)./sum(iX~=(Nstates+1),2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pauly Contracts 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PKH = PKH_prices_v4(CAPH,MU,rho,T);
Ppaid_pauly1 = PKH(iX(1,:),1)';
Ppaid_pauly2N = repmat(PKH(1,2:end)',1,nInd);
Ppaid_pauly = cat(1,Ppaid_pauly1,Ppaid_pauly2N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VIII. Welfare
%%%%%%%%%%%%%%%%%%%%%%%%%%%

  
%%%%%%%%%%%%%%%%%
% income %%%%%%%%
%%%%%%%%%%%%%%%%%

    age = [25:(25+T-1)]';
    
    file = dir(income_file);
        
    f = fullfile(file.folder,file.name);
    data=importdata(f);
    logincome=data.data;
    Wio_abi =exp(logincome(:,3)); 
    Wio_rea =exp(logincome(:,2)); 
    Wio_abi = Wio_abi(1:T,:);
    Wio_rea = Wio_rea(1:T,:);
       
    % flat income equal to NPV of real (representative)
    % note here there is no volatily
    D = mean(iX~=(Nstates+1),2).*power(rho,(1:T)-1)';
    NPVWio_rea = sum(D.*Wio_rea,1);
    Wio_rea_flat = ones(T,1)*NPVWio_rea/sum(D,1);
    
    % change dimensions, thousands of dollars
    Wio_rea_flat= Wio_rea_flat/1000;
    Wio_rea= Wio_rea/1000;
    Wio_abi = Wio_abi/1000;
    
    % first best 

%    Num_rea = mean(sum((Wio_rea_flat-repmat(Ppaid_spot_mean,1,nInd)).*(iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    Num_rea = mean(sum((Wio_rea-repmat(Ppaid_spot_mean,1,nInd)).*(iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    Den = mean(sum((iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    CFB_rea = Num_rea/Den*1000;

    Num_abi = mean(sum((Wio_abi-repmat(Ppaid_spot_mean,1,nInd)).*(iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    CFB_abi = Num_abi/Den*1000;

    % note, everything here in thousands of dollars
        
    [CE_N_abi,CE_ST_abi,CE_PKH_abi,CE_NBNS_abi,CE_HHW_abi,Cg_abi] = analysis_main(Wio_abi,CAPH,MU,Ppaid_N,Ppaid_spot,...
        Ppaid_pauly,Ppaid_spot_mean,iX,Nstates,rho,gamma,nInd,T,tolFP,maxreps);
    [CE_N_rea,CE_ST_rea,CE_PKH_rea,CE_NBNS_rea,CE_HHW_rea,Cg_rea] = analysis_main(Wio_rea,CAPH,MU,Ppaid_N,Ppaid_spot,...
        Ppaid_pauly,Ppaid_spot_mean,iX,Nstates,rho,gamma,nInd,T,tolFP,maxreps);   
    
    maintable = array2table([CFB_rea CE_ST_rea CE_HHW_rea CE_N_rea ;...
    CFB_abi CE_ST_abi CE_HHW_abi CE_N_abi]); 
    maintable = renamevars(maintable,["Var1","Var2","Var3","Var4"],...
                        ["CE_FB","CE_ST","CE_GHHW","CE_GLTHI"]);

    maintable.gainsGHHW = (maintable.CE_GHHW-maintable.CE_ST)./(maintable.CE_FB-maintable.CE_ST);
    maintable.gainsGLTHI = (maintable.CE_GLTHI-maintable.CE_ST)./(maintable.CE_FB-maintable.CE_ST);
    maintable.metric1 = maintable.gainsGLTHI./maintable.gainsGHHW;
    maintable.metric2 = (maintable.CE_GHHW-maintable.CE_GLTHI)./(maintable.CE_GHHW);

    simulation_output = struct();
    simulation_output.maintable = maintable;
    simulation_output.CE_PHI_abi = CE_N_abi;
    simulation_output.CE_PHI_rea = CE_N_rea;
    simulation_output.CE_HHW_abi = CE_HHW_abi;
    simulation_output.CE_HHW_rea = CE_HHW_rea;
    simulation_output.CE_PKH_abi = CE_PKH_abi;
    simulation_output.CE_PKH_rea = CE_PKH_rea; 
    simulation_output.Pspot_N = Pspot_N;
    simulation_output.Ppaid_spot = Ppaid_spot;
    simulation_output.Cg_abi = Cg_abi;
    simulation_output.Cg_rea = Cg_rea;
    simulation_output.iX = iX;
    simulation_output.Wio_abi = Wio_abi;
    simulation_output.Wio_rea = Wio_rea;
    simulation_output.Ppaid_N = Ppaid_N;
    simulation_output.Ppaid_pauly = Ppaid_pauly;
    simulation_output.CE_NBNS_abi = CE_NBNS_abi;
    simulation_output.CE_NBNS_rea = CE_NBNS_rea;
    simulation_output.PKH = PKH; 
    
    save(model_parameters.outputfile,'simulation_output');
         
end
        
    
